﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace WindGoes6.NumCS
{
	/// <summary>
	/// 矩阵类。
	/// </summary>
	public class Matrix
	{
		/// <summary>
		/// 矩阵行数，未初始化时为-1。
		/// </summary>
		public int RowCount { get; private set; } = -1;

		/// <summary>
		/// 矩阵列数，未初始化前为-1。
		/// </summary>
		public int ColumnCount { get; private set; } = -1;


		double[,] mdata { get; set; }


		/// <summary>
		/// 使用矩阵尺寸初始化矩阵。
		/// </summary>
		/// <param name="rowCount">行数。</param>
		/// <param name="columnCount">列数。</param>
		public Matrix(int rowCount, int columnCount)
		{
			this.RowCount = rowCount;
			this.ColumnCount = columnCount;
			mdata = new double[rowCount, columnCount];
		}

		/// <summary>
		/// 对输入矩阵进行复制，生成新矩阵。
		/// </summary>
		/// <param name="m">待复制矩阵。</param>
		public Matrix(Matrix m)
		{
			this.RowCount = m.RowCount;
			this.ColumnCount = m.ColumnCount;
			mdata = new double[m.RowCount, m.ColumnCount];
			for (int i = 0; i < m.RowCount; i++)
				for (int j = 0; j < m.ColumnCount; j++)
					mdata[i, j] = m.GetValue(i, j);
		}

		/// <summary>
		/// 对矩阵进行逐单元处理的函数。
		/// </summary>
		/// <param name="f">处理函数。</param>
		public void Process(Func<double, double> f)
		{
			for (int row = 0; row < RowCount; row++)
			{
				for (int col = 0; col < ColumnCount; col++)
				{
					SetValue(row, col, f(GetValue(row, col)));
				}
			}
		}


		/// <summary>
		/// 获得指定行列的值。注：为效率不做边界检测。
		/// </summary>
		/// <param name="row">指定的行。</param>
		/// <param name="col">指定的列。</param>
		/// <returns></returns>
		public double GetValue(int row, int col)
		{
			return mdata[row, col];
		}

		/// <summary>
		/// 输入相应的值，对矩阵进行设置。注：为了效率，不对边界进行检测。
		/// </summary>
		/// <param name="row">指定的行。</param>
		/// <param name="col">指定的列。</param>
		/// <param name="v">输入的值。</param>
		public void SetValue(int row, int col, double v)
		{
			mdata[row, col] = v;
		}

		/// <summary>
		/// 矩阵与常量加法计算。
		/// </summary>
		/// <param name="matrix">输入矩阵</param>
		/// <param name="operand">增量</param>
		/// <returns></returns>
		public static Matrix operator +(Matrix matrix, double operand)
		{
			Matrix m = new Matrix(matrix);
			for (int row = 0; row < m.RowCount; row++)
			{
				for (int col = 0; col < m.ColumnCount; col++)
				{
					m.SetValue(row, col, m.GetValue(row, col) + operand);
				}
			}
			return m;
		}

		/// <summary>
		/// 矩阵与常量减法计算。
		/// </summary>
		/// <param name="matrix">输入矩阵</param>
		/// <param name="operand">增量</param>
		/// <returns></returns>
		public static Matrix operator -(Matrix matrix, double operand)
		{
			Matrix m = new Matrix(matrix);
			for (int row = 0; row < m.RowCount; row++)
			{
				for (int col = 0; col < m.ColumnCount; col++)
				{
					m.SetValue(row, col, m.GetValue(row, col) - operand);
				}
			}
			return m;
		}

		/// <summary>
		/// 向控制台格式化输出矩阵。
		/// </summary>
		public void Print()
		{
			for (int i = 0; i < RowCount; i++)
			{
				for (int j = 0; j < ColumnCount; j++)
				{
					Console.Write("\t" + mdata[i, j]);
				}
				Console.Write("\n");
			}
		}


		/// <summary>
		/// 矩阵转置。
		/// </summary>
		/// <param name="Ma"></param>
		/// <returns></returns>
		public Matrix T()
		{
			Matrix mt = new Matrix(ColumnCount, RowCount);
			for (int i = 0; i < ColumnCount; i++)
				for (int j = 0; j < RowCount; j++)
					mt.mdata[i, j] = mdata[j, i];
			return mt;
		}

		//矩阵加法
		public static Matrix MatrixAdd(Matrix Ma, Matrix Mb)
		{
			int m1 = Ma.RowCount;
			int n1 = Ma.ColumnCount;
			int m2 = Mb.RowCount;
			int n2 = Mb.ColumnCount;
			if ((Ma.RowCount != Mb.RowCount) || (Ma.ColumnCount != Mb.ColumnCount))
			{
				Exception myException = new Exception("数组维数不匹配");
				throw myException;
			}

			Matrix Mc = new Matrix(Ma.RowCount, Ma.ColumnCount);
			double[,] c = Mc.mdata;
			double[,] a = Ma.mdata;
			double[,] b = Mb.mdata;

			for (int i = 0; i < Ma.RowCount; i++)
				for (int j = 0; j < Ma.ColumnCount; j++)
				{
					double temp = Ma.GetValue(i, j) + Mb.GetValue(i, j);
					Mc.SetValue(i, j, temp);
				}
			return Mc;
		}


		//矩阵减法
		public static Matrix MatrixSub(Matrix Ma, Matrix Mb)
		{
			if ((Ma.RowCount != Mb.RowCount) || (Ma.ColumnCount != Mb.ColumnCount))
			{
				Exception myException = new Exception("数组维数不匹配");
				throw myException;
			}

			Matrix Mc = new Matrix(Ma.RowCount, Ma.ColumnCount);
			double[,] c = Mc.mdata;
			double[,] a = Ma.mdata;
			double[,] b = Mb.mdata;

			for (int i = 0; i < Ma.RowCount; i++)
				for (int j = 0; j < Ma.ColumnCount; j++)
				{
					double temp = Ma.GetValue(i, j) - Mb.GetValue(i, j);
					Mc.SetValue(i, j, temp);
				}
			return Mc;
		}
		//矩阵乘法
		public static Matrix MatrixMulti(Matrix Ma, Matrix Mb)
		{
			int m1 = Ma.RowCount;
			int n1 = Ma.ColumnCount;
			int m2 = Mb.RowCount;
			int n2 = Mb.ColumnCount;
			if (Ma.ColumnCount != Mb.RowCount)
			{
				Exception myException = new Exception("数组维数不匹配");
				throw myException;
			}

			Matrix Mc = new Matrix(Ma.RowCount, Mb.ColumnCount);
			double[,] c = Mc.mdata;
			double[,] a = Ma.mdata;
			double[,] b = Mb.mdata;
			for (int i = 0; i < Ma.RowCount; i++)
				for (int j = 0; j < Mb.ColumnCount; j++)
				{
					c[i, j] = 0;
					for (int k = 0; k < n1; k++)
					{
						c[i, j] += a[i, k] * b[k, j];
					}
					//double temp = Ma.GetNum(i, j) * Mb.GetNum(j, i);
					// Mc.SetNum(i, j, temp);
				}
			return Mc;
		}
		//矩阵数乘
		public static Matrix MatrixSimpleMulti(double k, Matrix Ma)
		{
			Matrix Mc = new Matrix(Ma.RowCount, Ma.ColumnCount);
			//double[,] c = Mc.Detail;
			//double[,] a = Ma.Detail;

			for (int i = 0; i < Ma.RowCount; i++)
				for (int j = 0; j < Ma.ColumnCount; j++)
				{
					double temp = Ma.GetValue(i, j) * k;
					Mc.SetValue(i, j, temp);
				}
			return Mc;
		}


		//矩阵求逆（高斯法）
		public static Matrix MatrixInv(Matrix Ma)
		{
			int m = Ma.RowCount;
			int n = Ma.ColumnCount;
			if (m != n)
			{
				Exception myException = new Exception("数组维数不匹配");
				throw myException;
			}
			Matrix Mc = new Matrix(m, n);
			double[,] a0 = Ma.mdata;
			double[,] a = (double[,])a0.Clone();
			double[,] b = Mc.mdata;

			int i, j, row, k;
			double max, temp;

			//单位矩阵
			for (i = 0; i < n; i++)
			{
				b[i, i] = 1;
			}
			for (k = 0; k < n; k++)
			{
				max = 0; row = k;
				//找最大元，其所在行为row
				for (i = k; i < n; i++)
				{
					temp = Math.Abs(a[i, k]);
					if (max < temp)
					{
						max = temp;
						row = i;
					}
				}
				if (max == 0)
				{
					Exception myException = new Exception("没有逆矩阵");
					throw myException;
				}
				//交换k与row行
				if (row != k)
				{
					for (j = 0; j < n; j++)
					{
						temp = a[row, j];
						a[row, j] = a[k, j];
						a[k, j] = temp;

						temp = b[row, j];
						b[row, j] = b[k, j];
						b[k, j] = temp;
					}
				}
				//首元化为1
				for (j = k + 1; j < n; j++)
					a[k, j] /= a[k, k];
				for (j = 0; j < n; j++)
					b[k, j] /= a[k, k];

				a[k, k] = 1;
				//k列化为0
				//对a
				for (j = k + 1; j < n; j++)
				{
					for (i = 0; i < k; i++)
						a[i, j] -= a[i, k] * a[k, j];
					for (i = k + 1; i < n; i++)
						a[i, j] -= a[i, k] * a[k, j];
				}
				for (j = 0; j < n; j++)
				{
					for (i = 0; i < k; i++)
						b[i, j] -= a[i, k] * b[k, j];
					for (i = k + 1; i < n; i++)
						b[i, j] -= a[i, k] * b[k, j];
				}
				for (i = 0; i < n; i++)
					a[i, k] = 0;
				a[k, k] = 1;
			}
			return Mc;
		}
		//矩阵求逆（伴随矩阵法）
		public static Matrix MatrixInvByCom(Matrix Ma)
		{
			double d = Matrix.MatrixDet(Ma);
			if (d == 0)
			{
				Exception myException = new Exception("没有逆矩阵");
				throw myException;
			}
			Matrix Ax = Ma.AdjointMatrix();
			Matrix An = Matrix.MatrixSimpleMulti((1.0 / d), Ax);
			return An;
		}

		//对应行列式的代数余子式矩阵
		public static Matrix MatrixSpa(Matrix Ma, int ai, int aj)
		{
			int m = Ma.RowCount;
			int n = Ma.ColumnCount;
			if (m != n)
			{
				Exception myException = new Exception("数组维数不匹配");
				throw myException;
			}
			int n2 = n - 1;
			Matrix Mc = new Matrix(n2, n2);
			double[,] a = Ma.mdata;
			double[,] b = Mc.mdata;

			//左上
			for (int i = 0; i < ai; i++)
				for (int j = 0; j < aj; j++)
				{
					b[i, j] = a[i, j];
				}
			//右下
			for (int i = ai; i < n2; i++)
				for (int j = aj; j < n2; j++)
				{
					b[i, j] = a[i + 1, j + 1];
				}
			//右上
			for (int i = 0; i < ai; i++)
				for (int j = aj; j < n2; j++)
				{
					b[i, j] = a[i, j + 1];
				}
			//左下
			for (int i = ai; i < n2; i++)
				for (int j = 0; j < aj; j++)
				{
					b[i, j] = a[i + 1, j];
				}
			//符号位
			if ((ai + aj) % 2 != 0)
			{
				for (int i = 0; i < n2; i++)
					b[i, 0] = -b[i, 0];
			}
			return Mc;
		}


		//矩阵的行列式
		public double Det(Matrix Ma)
		{
			int m = Ma.RowCount;
			int n = Ma.ColumnCount;
			if (m != n)
			{
				Exception myException = new Exception("数组维数不匹配");
				throw myException;
			}
			double[,] a = Ma.mdata;
			if (n == 1)
				return a[0, 0];
			double D = 0;
			for (int i = 0; i < n; i++)
			{
				D += a[1, i] * MatrixDet(MatrixSpa(Ma, 1, i));
			}
			return D;
		}

		//矩阵的行列式
		public static double MatrixDet(Matrix Ma)
		{
			int m = Ma.RowCount;
			int n = Ma.ColumnCount;
			if (m != n)
			{
				Exception myException = new Exception("数组维数不匹配");
				throw myException;
			}
			double[,] a = Ma.mdata;
			if (n == 1)
				return a[0, 0];
			double D = 0;
			for (int i = 0; i < n; i++)
			{
				D += a[1, i] * MatrixDet(MatrixSpa(Ma, 1, i));
			}
			return D;
		}

		/// <summary>
		/// 返回当前矩阵的伴随矩阵(adjoint matrix).
		/// </summary>
		/// <returns></returns>
		public Matrix AdjointMatrix()
		{
			Matrix m = new Matrix(RowCount, ColumnCount);
			for (int row = 0; row < RowCount; row++)
				for (int col = 0; col < ColumnCount; col++)
					m.mdata[row, col] = MatrixDet(MatrixSpa(this, col, row));
			return m;
		}

		/// <summary>
		/// 两个矩阵相加。
		/// </summary>
		/// <param name="ma">第一个加数</param>
		/// <param name="mb">第二个加数</param>
		/// <returns></returns>
		public static Matrix operator +(Matrix ma, Matrix mb)
		{
			return MatrixAdd(ma, mb);
		}

		//重载减号
		public static Matrix operator -(Matrix Ma, Matrix Mb)
		{
			return MatrixSub(Ma, Mb);
		}

		/// <summary>
		/// 以二维数据的形式输出矩阵
		/// </summary>
		/// <returns></returns>
		public override string ToString()
		{
			StringBuilder sb = new StringBuilder();
			for (int i = 0; i < RowCount; i++)
			{
				for (int j = 0; j < ColumnCount; j++)
				{
					sb.Append(mdata[i, j].ToString() + "\t");
				}
				sb.AppendLine();
			}
			return sb.ToString();
		}

	}
}
